Dynamics of the inhomogeneous Dicke model 
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We study the time dynamics of a single boson coupled to a bath of two-level systems (spins 
1/2) with different excitation energies, described by an inhomogeneous Dicke model. Analyzing 
the time-dependent Schrodinger equation exactly we find that at resonance the boson decays in 
time to an oscillatory state with a finite amplitude characterized by a single Rabi frequency if 
the inhomogeneity is below a certain threshold. In the limit of small inhomogeneity, the decay is 
suppressed and exhibits a complex (mainly Gaussian-like) behavior, whereas the decay is complete 
and of exponential form in the opposite limit. For intermediate inhomogeneity, the boson decay is 
partial and governed by a combination of exponential and power laws. 

PACS numbers: 42.50.-p,71.35.Lk,06.20.-f,03.67.-a 
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I, INTRODUCTION 

Coherent interaction between light and matter [lj con- 
tinues to receive strong interest due to significant exper- 
imental progress in various areas of physics. Prime ex- 
amples are the achievement of Bose-Einstein condensa- 
tion of cold-atom gases in electromagnetic traps [2] which 
made possible the coherent coupling of 10 5 atoms to a 
single photon of an optical resonator 0, 0]- The time 
dynamics of quantum optical systems has received par- 
ticular attention @, [f| due to fast optical probing tech- 
niques, especially in the context of quantum metrology 
based on cavity-QED systems containing atomic ensem- 
bles @, @|- Advances in solid-state technology enabled 
the fabrication of optical microcavities in semiconduc- 
tors where electron-hole excitations in quantum wells are 
str ong ly coupled to a photon eigenmode of the cavity 
[3, llOj . Strong coupling of a transmission- line resonator 
to a Cooper-pair box as well as coupling of a cav- 
ity to a single semiconductor quantum dot have been 
demonstratedfl2l. [l3| . Several schemes for quantum com- 
puting based on light-matter interaction have been pro- 
posed [13, [H, m, E3, [ili- 

The theoretical understanding of all these coupled 
light-matter systems is based on a model introduced long 
ago by Dicke pjj], which describes N two-level systems 
('spin bath') with excitation energies ej coupled to a sin- 
gle boson mode u> of the quantized light-field, see Eq. ([I]) 
below. For the special case of identical atoms (ej = e) 
and constant couplings constant gj between boson and 
spin bath this model has been diagonalized [2fj| . and the 
time dynamics obtained exactly [21(. For inhomogeneous 
gj (but still constant ej) the boson was shown to oscil- 
late with a single Rabi frequency £1 = y/N (g 2 ), where 
\/ (g 2 ) is an effective spin coupling. Also perturbative 
[24] and numerical [23| approaches to the time dynamics 
were considered. 

In this paper we solve the quantum time-dynamics of 
a single boson mode coupled to a bath of non-identical 
spins 1/2 characterized by inhomogeneous energy ('Zee- 
man') splittings ej with bandwidth A. In condensed mat- 



ter systems such energy inhomogeneities are generally ex- 
pected, a typical example being the exciton-polariton sys- 
tem where such inhomogeneities arise from the unavoid- 
able disorder in a semiconductor [HI]. In quantum optical 
systems atomic levels are usually quite perfect (ej = e); 
however, for example, in cold-atom QED systems such 
inhomogeneities can play a role as a trap induces spatial 
variation of the magnetic field [HI] . 

Analyzing the time-dependent Schrodinger equation 
exactly we find that the bosonic occupation number de- 
cays only partially if the inhomogeneity A is below a 
threshold given by a single Rabi frequency fi. Below the 
threshold the boson decays to an oscillatory state deter- 
mined by il and a reduced amplitude which decreases 
with increasing ratio A/Q. The time decay is exponen- 
tial for large spin-bath inhomogeneity A ^> O, is complex 
(mainly Gaussian-like) in the opposite limit A <C £1 and 
is a combination of exponential and power law behavior 
in the intermediate regime A ~ f2. These results are 
valid if the boson energy is tuned in resonance with the 
average spin excitation energy (e) — u> = 0. With increas- 
ing detuning |(e) — u)\ ^> max{il, A} the time dynamics 
of the boson becomes suppressed. 

The paper is organized as follows. In Section[Ulwe ana- 
lyze the time-dependent Schrodinger equation and derive 
the exact solution in the Laplace domain. In Section Hill 
we consider rectangular and Gaussian distribution func- 
tions of €j in resonance with the boson mode, (e) = uj, 
to obtain the time evolution of the wave functions. Sec- 
tion IIVI contains the analysis and discussion of a finite 
detuning, (e) ^ id. In the appendices we give details of 
the calculations in Sections [Ul] and [ 



II, THE INHOMOGENEOUS DICKE MODEL 

The Hamiltonian for the Dicke model governing the 
dynamics of a single boson mode coupled to N two-level 
systems is given by 



JV 



N 



H = cotfb + £ ej S? + 9 j {S+b + S-tf) , (1) 



3=1 



2 



5? ± iS?, and 



where are spin-1/2 operators, 

b (tf) the standard Bose annihilation (creation) operator 
|27l |. The total number of excitations, L = n + J2j S*, 
is conserved in the Dicke model, where n = b'b is the 
bosonic occupation number. The eigenvalues c of L are 
the so-called cooperation numbers, given by c = (L), 
where (...) denotes the expectation value. 

In the following we assume that the spin bath can be 
prepared in its ground state with all spins down, e.g. 
either dynamically or by thermal cooling [28]. Also, the 
mode u is assumed to be empty or occupied by one boson 
only. The non-equilibrium dynamics of a single boson ex- 
citation can then be initiated by a short radiation pulse 
from an external source. The dissipation of the boson 
mode, e.g. through leakage of photons through the mir- 
rors that define an optical cavity, can be used to detect 
the dynamics if the cavity escape time exceeds the inter- 
nal time scales of the system dynamics. In a multi-shot 
experiment @, @] the probability of detecting a leaking 
photon at a given time is proportional to the boson ex- 
pectation value. Next, we note that if initially the system 
has only one excitation, either in the spin or in the boson 
subsystem, the subsequent time evolution is restricted to 
this subspace and described by the general state 



\9(t))=a(t)\^,l) + ^20 j (t)\^j,O) 



(2) 



3=1 



with c = —N/2+1, and where a (t) and (3j (t) are normal- 
ized amplitudes, \a (t)\ 2 + J^j \Pj\ 2 — lj °f finding either 
a state with one boson and no spin excitations present or 
a state with no boson and the j th -spin excited (flipped) 

The time evolution within this subspace is determined 
by the interaction term in Eq. (HJ that transfers back 
and forth the excitations between the spin bath and 
the boson. Inserting |\I> (i)) into the time-dependent 
Schrodinger equation we obtain 



da (t) 

- 1 — ; — 
dt 

d(3 k (t) 

i ; 

dt 



E 



(3) 



In above derivation we have subtracted the integral of 
motion u>L from the Hamiltonian Eq. ([I]) as it leads only 
to an overall phase of \\V) with no observable effect. The 
initial conditions a (0) = 1, /3j (0) = assumed in the 
following correspond to a singly occupied boson mode. 
The physical observable of interest is the time-dependent 
expectation value of the boson occupation number, which 
can be expressed in terms of the amplitude a as (n (t)) — 

(*(t)|n|* (*)> = !«(*) I 2 - 

The set of equations, Eq. ((SJ, is equivalent to the one 
obtained in the Weisskopf-Wigner theory in the study 
of bosonic systems in contrast of spins 1/2 con- 
sidered here. We solve Eq. ((3j) by making use of the 



Laplace transform, a (s) = f °° dta (t) e~ st , 3fJs > 0. In 
the Laplace domain we obtain then a system of linear 
algebraic equations. By solving them we find 



a(s) 



s+(w— e.,-)iV/2— w+£j 



(4) 



where (...) = ■ ■ - )/N, The sum over j depends 
on the particular form of the inhomogeneities of tj and 
cjj. To be specific, we consider the following limiting 
cases when 6j varies on a much longer or shorter length 
scale than gj, which also includes the case with either €j 
or gj being constant. In this case and for large N the 
sum can be substituted by an integral, ■ ■ ■ )/N — * 
J dedg P(e)Q(g), where P (e) and Q (g) are indepen- 
dent normalized distribution functions of the excitation 
energies and coupling constants, respectively. The inte- 
gral over g in Eq. (@1) separates and gives an effective 
coupling \J (g 2 ) [31| . Further, we assume that the boson 
mode u) is tuned in resonance with the spin bath, i.e. 
u-(e) = 0. 



III. INVERSE LAPLACE TRANSFORM 

The inverse Laplace transform of Eq. ([4]) depends on 
the particular form of P (e) that determines the ana- 
lytic structure of a(s). We will analyze several cases 
below. If the spin bath is homogeneous then P (e) = 
S (e — uj), and a (s) has two poles on the imaginary axis 
at s = ±iyj N (g 2 ), with the associated residues 1/2. In 
the time domain these poles give a (t) = cos(ftt), where 
^ = y/N (d 2 ) is the collective Rabi frequency due to all 
N spins. This agrees with the result obtained from exact 
diagonalization [2l| . 

Next, we consider an inhomogeneous spin bath with ex- 
citation energies spread over a band of width A, for which 
we have P (e) = 9 (— e + ui + A/2) 8 (e - u> + A/2) /A, 
where 9 (x) is the step function. This case is realized 
e.g. for e 3 = jA/N, -N/2 < j < N/2, i.e. spins in a 
magnetic field with constant gradient. The integral over 
e in Eq. Q gives 



a(s) 



A 



In 



/ is-A/2 \ 



(5) 



Note that the inverse Laplace transform of Eq. Q is 
in principle a quasi-periodic function of t. Therefore, Eq. 
© is correct up to the Poincare recurrence time t p which 
we can estimate as follows. We evaluate the discrete sum 
over 6j exactly, expand it in 1/N, and estimate the time 
at which corrections to the logarithmic term in Eq. © 
(due to discretness of the sum) become important to be 
t p = N/A. Thus, the following time behavior is valid 
for times less than t p = N/A. For small N it is more 
convenient to find the few poles of Eq. (@]) directly and 
analyze a(t) numerically as a sum of few harmonic modes 
rather than to use Eq. ([5]). 



We discuss now the analytic structure of a (s) in Eq. 
©. There are two branch points at s = ±iA/2 due to 
the logarithm. We choose the branch cut as a straight 
line between these two points. In addition, there are two 
poles at s — ±zso given by the zeroes of the denominator 
where sq is a real and positive solution of 



exp 



s A 



so - A/2 
s + A/2' 



(6) 



In the time domain, the amplitude a has two contribu- 
tions, a — a p + a c . One is given by the poles, 

a * {t) = l + iV(.g 2 )/ 2 (^-AV4) C ° S M) ■ (?) 

This contribution describes a residual oscillation at long 
times with amplitude that is reduced from the initial 
value a (0) = 1. The other one is given by the integral 
enclosing the branch cut, 



a c (t) 



dy 



(W(g 2 )/A 2 ) cos(2/At/2) 



y 



A 5 "" 



In (SO 



/V<5 2 



(8) 

This contribution describes the decay that occurs due to 
destructive interference of many modes forming a contin- 
uous spectrum (for large N). 

The integral in Eq. ([5]) can be approximated quite 
accurately for t ^ 2/ A. Due to the fast oscillating cosine 
the main contribution to the integral comes from y < 
2/ At <C 1. Expansion of the logarithm in Eq. (JSj) for 
small y permits us to evaluate the integral in terms of the 
Integral Sine and Cosine. An expansion of these special 
functions for At/ 2 ^ 1 gives 



•.(*) = 



A 2 



N(g 2 



Ae -AAt/2 

2^ 



A 2 sin 



(At/2) \ 



7T 2 (1+A 2 ) At/2 J 



(9) 



where A = tt/2/ |l - A 2 /AN (g 2 )\. Note that for vanish- 
ing coupling g, a p (t) vanishes and a c (t) tends to one. 
Further, the integrand in Eq. ([5]) can be expanded for 
A 2 /N (g 2 ) for A 2 < N (g 2 ). The leading term is lin- 
ear in A 2 /N (g 2 ) and the remaing integral in the pref- 
actor is a complicated decaying function of t which we 
approximate qualitatively. First, we perform a change 
of variable - y = tanh(x) turning the denominator into 
1/ / = cxp (— log (/)), where we expand log (/) up to x 2 
and linearize the argument of the cosine in x for x <C 1 . 
Finally, as a result of the Gaussian integral over x we 
obtain a Gaussian decay law, 



a c (i) 



2nN (g 2 ) V 7T 2 + 4 



■ exp 



7T 2 A 2 t 2 

4(tt 2 +4) 



(10) 



This approximation agrees reasonably well with Eq. JH]) 
when evaluated numerically for t < 6/ A but breaks down 
for t > 6/A where Eq. ([9]) is valid, see Fig. [Q 

The time-dynamics of (n (f)) = |a(i)| 2 can be clas- 
sified in terms of the ratio f2/A, with Rabi frequency 




A?r/ Sq 



At/2 



Figure 1: Time evolution of the boson (n(fj) — |o:(t) | 2 ob- 
tained from numerical evaluation of Eqs. (JXl E|> - full lines. 
Period of oscillation is T — 2n/so and grey bars are |a p (0)| 2 . 
Main plot: A/fi = 2.2; dashed line: \a p (t) + a c (t)\ 2 fro m 
Eqs. Q i, s o = 0.57A, and |a P (0)| 2 = 0.26. Inset: 
A/A = 0.2; dashed line: \a p (0) ± a c (t)\ 2 from Eqs. iffl [Tojl 
(valid for At/2 < 3) and from Eq. © (valid for At/ 2 > 3). 
The decay of |a(i)| 2 is small, (|a P (0)| 2 = 1 - A 2 /6fi 2 ), and 
the main contribution comes from a„. 



^ = y/ N (g 2 ) . If the inhomogeneity of the spin bath 
is small, A fi, the boson oscillates with a single fre- 
quency like in the homogeneous case. The main con- 
tribution to a (t) comes from the poles Eq. ([7]) with 
so = fi+ A 2 /24fi, which is shifted with respect to the ho- 
mogeneous system. The amplitude of a (t) is only slightly 
reduced from its initial value, 1 — A 2 /12f2 2 . The decay 
law to this value is mainly Gaussian-like, Eq. (flOl) . with 
the decay time t\ ~ 2. 4/ A, see Fig. [TJ If the spin bath 
is strongly inhomogeneous, A 3> ft, the boson mode de- 
cays completely from a (0) = 1 to 0. The main contri- 
bution to a (t) comes from the branch cut, Eq. ([9]), with 
A w 27r£! 2 /A 2 , whereas the pole contribution is exponen- 
tially small. The decay behavior is mainly exponential 
with timescale £2 ~ A/ttQ 2 . At long times t ^> ti the 
second term in Eq. |(9]) becomes dominant, exhibiting a 
slow power-law decay. 

In the intermediate regime, A ~ f2, the time decay is 
only partial, with the amplitude of the residual oscillation 
of a (t) being less than unity but staying constant in time. 
Its precise value can be found from the numerical solution 
of Eqs. (|6|7p . The decay displayed in Eq. ((9]) is governed 
by a combination of exponential and power law behavior. 
As A ~ 1 and so — & — A there is no clear separation 
of time scales coming from the exponential, the inverse 
power law and the oscillatory contribution, see Fig. [TJ 
Note that in case of A = 2f2 the first term in Eq. © 
vanishes, thus the decay in this particular case is purely 
power law. The non-standard dynamics, in particular 



4 



the non-exponetial decay in the intermediate regime, is a 
manifestation of the quantum nature of the system. For 
other models with non-Markovian decay see e.g. [3cLl3^|. 

For a Gaussian distribution P (e) = 
exp(— (e — ui) 2 /A 2 )/y / 7tA the dynamics we find is 
qualitatively the same as the one obtained before for the 
rectangular distribution, see Appendix A. The e-integral 
in Eq. (@]) leads to the complex error function of s. 
In Laplace space, a (s) exhibits one branch cut along 
the imaginary axis that vanishes at ±zoo. In the time 
domain, a (t) is given by an integral around this branch 
cut. In the limit of A <C we recover the previous 
result for the homogeneous spin bath. In the opposite 
limit of strong inhomogeneity, A ^> ft, we obtain the 
same result as in Eq. ^ up to numerical prefactors y/ir. 

The physical interpretation of the decay is as follows. 
The boson flips, say, spin j, and then this spin precesses 
for some time with a frequency ej before this excitation 
gets transfered back to the boson. The acquired phase 
of the boson is thus different for each particular spin. 
The sum over these random phases eventually leads to 
destructive interference (for N ^> 1) and thus to a decay. 



coupling can be removed to lowest order in g and thereby 
an effective XY spin-coupling within the spin bath is ob- 
tained (lq . As a result, the boson number n and the 
z-component of the total spin J2j $j are conserved sep- 
arately by this effective Hamiltonian. Thus, again, the 
initially excited boson mode will remain unaltered under 
the evolution in leading order of the perturbation. How- 
ever, there is a virtual boson process which induces the 
dynamics within the spin bath. 

V. CONCLUSIONS 

In conclusion, we analyzed the dynamics of a single bo- 
son mode coupled to an inhomogeneous spin bath exactly, 
and found a complex decay behavior of the boson. While 
we focused in this work on particular inhomogeneities of 
the spin bath excitation energies, it is straightforward to 
apply the approach presented here to other cases. 
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IV. FINITE DETUNING 

Next, we analyze the effect of finite detuning. If the 
spin bath is homogeneous, a small detuning \oj — e| <C 
ft forces a (t) to oscillate with two distinct frequencies 



(N - 1) (w - e) /2± \J (e - uj) + instead of only one 
ft. A large detuning \u> — e\ » ft suppreses the dynamics 
of a (t) . The phase of the wave function oscillates with 
frequency N (u> — e) /2 but the amplitude stays constant 
at the initial value of a (0) = 1 up to a small correction 
of the order of ft 2 / (ui — e) 2 . 

In the inhomogeneous case we perform a similar cal- 
culation as for zero-detuning and obtain a (s) with an 
analytic structure similar to Eq. |j5j, see Appendix B. 
There are two poles on the imaginary axis and a branch 
cut that is responsible for the relaxation. Explicit ex- 
pressions for a p (t) and a c (t) can be obtained and are 
generalizations of Eqs. (|7|8|l . see Eq. IjBip . For small 
detuning |(e) — u\ -C ft two poles emerge that are not 
complex conjugates of each other and thus lead to two 
distinct frequencies of the residual oscillations of a (t). 
Large detuning |(e) — u>\ 3> max{ A, ft} suppresses the 
relaxation and any long-time dynamics. The main con- 
tribution to a (t) comes from one of the poles with residue 
1 - ft 2 / (u - (e)) 2 . Thus, the initial value a (0) = 1 re- 
mains almost unaltered under evolution independent of 
the ratio ft/ A. 

The dynamics at large detuning can also be analyzed 
using perturbation theory. Applying a Schrieffer- Wolff 
transformation to the Dicke Hamiltonian the boson-spin 
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Appendix A: GAUSSIAN DISTRIBUTION OF 
THE SPINS' SPLITTING ENERGIES 



Here we derive the time dynamics resulting from 
the Gaussian distribution function of e, P (e) = 
_^_ e -(c-w) /A _ Performing the integral over e in Eq. 
(j4]) we obtain 



a(s) 



(Al) 



where zu (z) is defined in the upper and lower complex 
half-planes separately as 

, . I w (z) , Im z > ... 
™{z) = \ \! ' " (A2) 
I — w \—z) , Im z < (J 

and where w (z) — e~ 7 - 2 enic (— iz) is the error function. 
The function w (z) has a branch cut along the real axis, 
lirn^o m (±i<5) = ±1, which vanishes at infinity, 

2 

lim lim w (x ± i5) = lim e~ x (±1 + erf (ix)) = 0. 

x — »±oo (5 — >0 x — >±oo 

(A3) 

The inverse Laplace transform is given by an integral 
around the entire imaginary axis 
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N(o 2 ) r°° -iyt/2A -4y 2 

a (t) = -p-L / dy J- (A4) 

2 ^ A2J ~°° Ly-^£± W {2y)) + VE»m^ ( ly -^l w{ 2y)) 



where the substitutions s = t2Ay and uj (—z) = 2e 2 — 
u> (z) were used. 

For A < y/N {g 2 ) Eq. ((Alj) can be expanded in the 
small parameter A. The leading term has an analytical 
structure similar to Eq. |(5|). There are two symmetric 
poles on the imaginary axis and a finite length branch 
cut between s = ±z2A. The contribution from the poles 
is 



ot p (t) 



2 cos (sot) 



(A5) 



1-N{g*) /si 
Using the large z asymptotics of the error function 

_ 2 2 

erfc(z) = ^7(1 -222)1 the two poles are given by 



So = ±iy/N (g 2 ). The residues at this poles are 

e Sot 

Res a(s)e st = ——. 

s=s 2 



(A6) 



Thus, the contribution from the poles is dominant. In 
this limit we recover the non-interacting case, a single 
Rabi oscillation, 



:(t) =cos(0V( ff 2 } 



(A7) 



In the opposite regime A 3> \J N (g 2 ) there are no 
poles and there is just a single branch cut. The long 
time asymptotics can be evaluated by expanding the de- 
nominator for small y and approximating e~ Ay s» 1 in 
the numerator, 

N (9 2 ) f 1 , cos(yi/27) , k x 

«c(*) = ^=xr / dy ,J„.' 2 - (A8) 



V^A2 J 



This integral, up to a numerical factor, is the same as in 
Eq. flU in this limit. 



Appendix B: CALCULATION FOR <e> / uj 

Here we assume that the detuning is finite 7 = (e)— u> ^ 
0. We repeat the same steps as before and similarly to 
the zero detuning case we obtain in the Laplace domain 



a 00 = 



Nj_ , N(g2) -, / is+(N-2)~//2-A/2 
2 ^ A m Ws+(jV-2)7/2+A/2 



(Bl) 



This function is characterized by two poles and one 
branch cut. 

The two poles are given by zeroes of the denominator 
s = 1 (Nj/2 + si,2) where are the solutions of 



exp 



sA 



_ a -7- A/2 
N(g 2 ) ) ~~ S-7 + A/2' 



(B2) 



This equation is not symmetric with respect to s — > — s, 
thus the two poles are not symmetric. The residues of 
the poles are given by 



Resa s e s 



1 



N{a 2 ) 

(si,2-7) 2 -A 2 /4 



(B3) 



Performing the inverse Laplace transformation we obtain 
similarly to Eq. |(6|) 



•(*)= E 



e iNyt/2+is k t 

I ^(fl 2 ) ~ 

fe=l,2 1 ( Sfc - 7 ) 2 -AV4 



(B4) 



The branch points are s — 1 (JV7/2 ± A/2). Similarly 
to the case of zero detuning the contribution from the 
branch cut is given by the integral 



Ol c (t) 



2N{g 2 ) 



ol ((JV-2) 7 +A)t/2 



A 2 



dy 



iyAt/2 



tL 
A 



A s 



hi 



2irN(g 2 
A* 



(B5) 



At small detuning |w-(e)| < max (A/2, N (g 2 )) 
there are two distinct frequencies in Eq. (|F34[) . thus 
the final state oscillates with two frequencies. For a 



large detuning \u — (e)| » max (A/2, N (g 2 )) the relax- 
ation is suppressed. In the limit of strong detuning the 
roots of Eq. (|B2j) are given by Si = —2N(g 2 } and 
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S2 = 7 — A/2. The residue at S2 is exponentially small 
and the contribution from the poles is given only by the 
pole si 

-i2N(g 2 )t/A+iN-ft/2 
+ (2W( s 2 )/7+7 )^_ A 2/4 



In this result the amplitude of a (t) remains constant 
in time. From the initial condition a p (0) + a c (0) = 1 
the contribution from the branch cut is negligible, and 
therefore there is no decay for sufficiently strong detun- 
ing. The corrections to this result are small and of the 
order of max (A/2, N (g 2 )) / \u> - (e)|. 
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